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Abstract. We discuss methods for numerically solving the generalized Master equation GME which governs 
the time-evolution of the reduced density matrix p of a mechanically movable mesoscopic device in a dissipative 
environment. As a specific example, we consider the quantum shuttle - a generic quantum nanoelectromechanical 
system (NEMS). When expressed in the oscillator basis, the static limit of the GME becomes a large linear non- 
sparse matrix problem (characteristic size larger than 10 4 x 10 4 ) which however, as we show, can be treated using 
the Arnoldi iteration scheme. The numerical results are interpreted with the help of Wigner functions, and we 
compute the current and the noise in a few representative cases. 
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Introduction 



2. The generalized Master Equation (GME) 



Microelectromechanical systems (MEMS) are ap- 
proaching the nanoscale, which ultimately implies that 
the mechanical motion needs to be treated quantum 
mechanically. An example is the experiment by Park 
et al. 1 , where the break junction technique was 
used to create a single-electron transistor with a C60- 
molecule as the active part. The measured IV-curves 
display features that can be related to the mechani- 
cal vibrations of the molecule. Already earlier it was 
suggested theoretically 2 that a nanoscopic movable 
metallic grain, when operated in the Coulomb block- 
ade regime, can move charges one-by-one between 
source and drain contacts. This orderly transport was 
coined the shuttle regime. In recent years our group 
has developed theoretical methods to analyze the shut- 
tle transition in the quantum regime [31 El El E] , focus- 
ing not only on the IV-curve, but also considering noise 
and full counting statistics, which are important diag- 
nostic tools in unravelling the microscopic transport 
mechanisms. In this paper we present examples of our 
numerical results and discuss pertinent computational 



The Hamiltonian for the system includes terms de- 
scribing (i) the electronic part of the movable quantum 
dot (QD for short), (ii) its mechanical motion (which 
is quantized), (iii) the position dependent coupling of 
the QD and the leads, (iv) the leads (treated as nonin- 
ter acting fermions), and (v) coupling to environment, 
which damps the mechanical motion EI El E] • Using 
methods familiar from quantum optics, we integrate 
out the environmental degrees of freedom (the lead 
electrons, and a generic heat bath) to obtain a gen- 
eralized Master equation for the "system" (= QD + 
quantized oscillator) density operator: 



f)(t) = Cp(t) = (£ C oh + £driv + Aiamp)p(£). (1) 



Here £ co h,£driv an d £damp are superoperators corre- 
sponding to the coherent evolution, coupling to leads, 
and damping of the QD. It is sufficient to consider the 
diagonal electronic components (i.e., an empty and an 
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occupied QD, respectively), which satisfy j3J|Hl 

i r 2x 

poo(t) = t^[H osc , p 00 (t)] - —(e~~p 00 (t) 

+A)oW e_a? ) + T R e^pu(t)e^ 
amp Poo(^) , 

pn(t) = J^l H osc -eEx,pn(t)} + T L e~^ p 00 (t)e~ 
-^(e-p n (t)+pn(t)e-) 



+£damp Pll(t)- 



(2) 



where 



C damp p {P> P)] -^{N + 1/2) [x, [x, p)] . 

The physical parameters defining the quantum shut- 
tle are thus the couplings to leads T L / R , the oscillator 
frequency lj, the damping rate of the oscillator 7, the 
temperature T, and the tunnelling length A. 

3. The numerical solution 

We seek the steady state, poo = P11 = 0, i.e., the 
null vector of the Liouvillean, Cp = 0. Suppose we 
keep N lowest energy states of the oscillator. Since 
Poo/11 are full matrices in the oscillator basis, they 
both have TV 2 elements. The matrix representation of 
the Liouville superoperator has thus 27V 2 x 27V 2 ele- 
ments. (The situation is even more demanding if one 
considers a shuttle with more than just two electronic 
states, see 00.) We attempted the solution of these 
equations with standard Mat lab routines, such as the 
singular value decomposition, by gradually increasing 
the cut-off TV. As illustrated in Fig. d the Liouville 
matrices are not sparse. Unfortunately, reliable con- 
vergence could not be achieved once N approached few 
tens, which was quite inadequate for physical reasons, 
which indicated the necessity of using TV c± 100 • • • 200. 
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Figure 1: The nonzero matrix elements of the Liouville su- 
peroperator for a three QD device studied in |3 15] for three 
values of the cut-off N = 2,5, 10. nz is the number of non- 
zero elements. 



Clearly, more powerful numerical schemes are re- 
quired. One such method is the Arnoldi iteration [3J, 
which has clear advantages compared to singular value 
decomposition both in terms of computational speed 
and memory requirements. First, it is not necessary 
to store the full Liouvillean matrix, and, second, when 
one seeks the best approximation to the null vector it 
is possible to work in spaces which are much smaller 
than the full Liouville space. The central concept in 
the Arnoldi scheme is the Krylov space, defined as 



JCj = span(x , Lx , L J 1 x ), 



(3) 



where j is a small integer. The vector xo is a vec- 
torization of some arbitrary state represented by the 
two matrices Poo/11? i.e. a vector of length 27V 2 . Note 
that the calculation of the vectors I/xo, i = — 1) 
can be done directly from the GME, without storing 
the full Liouvillean matrix. The next step consists of 
finding an ort honor mal basis in the Krylov space, let 
us denote this by {qi},i = 1, j. The method pro- 
ceeds now by finding an approximate null vector in the 
Krylov space, spanned by the vectors {q^} thus involv- 
ing matrices of size j x j. In practice, we have found 
that j = 20 is sufficient, which makes the calculations 
very fast, and no supercomputing is necessary. Once 
the optimal vector in the Krylov space has been found, 
we can use the coordinates of this vector to construct 
the estimate xi (which is a vector of length 2N 2 ) from 
which the matrices Poo/11 can be reassembled. If the 
result is not sufficiently close to the null vector of the 
Liouvillean, one repeats the procedure by constructing 
the next Krylov space, now using xi as the seed. 

The Arnoldi scheme is iterative in its nature, and 
one must address the issue of convergence. For ex- 
ample, it is not a priori clear how many iterations are 
needed, and indeed convergence is not always achieved. 
This problem is solved by the use of a preconditioner. 
The basic idea is to find an invertible operator M 
in the Liouville space, such that the original problem 
£p stat = is cast in the form A^[£[p stat ]] = 0, and that 
the truncated version of the operator AiC gives rise 
to a rapidly converging iteration scheme. The Arnoldi 
scheme is particularly efficient in finding good approx- 
imations to eigenvalues (and corresponding eigenvec- 
tors) that are separated from the rest of the spectrum. 
Thus, in our case, the preconditioner should move the 
eigenvalues with a non- vanishing real part away from 
the origin of the complex plane. A good candidate for 
the preconditioner is to use the Sylvester part ^01 of 
C: 

£ [P] = £Sylv [p] + £rest [p] , (4) 



Modelling of Quantum Electromechanical Systems 3 



where the Sylvester part has the structure 

= ( -4ooPoo + Poo Ax) ^ 

where the elements A o/n can be gleaned off from the 
GME. The Sylvester part is rapidly invert ible, and we 
thus chose A4 = £gy lv . In general, it was essential 
to use the preconditioner, however we stress that the 
choice is somewhat subtle and there is no unique algo- 
rithm for this. We have encountered special situations 
where the Sylvester preconditioner did not work sat- 
isfactorily, and more work is required in refining the 
numerics in this case. 

Once the static density matrix is solved, the current 
is readily calculable from 

js tat = eTrosc {r R e 2 */ x pff} 

= eTV osc {r L e- 2 ^Voo at }- (5) 

We have recently shown that also the noise, and 
even the higher cumulants 0, can be calculated with 
similar methods. In particular, we find that the Fano 
factor F — S(0)/2eI (here S(0) is the zero- frequency 
component of the noise spectrum) can be expressed as 

F = 1 - [e* [OC-'C 

Here Q is a projection operator that projects away 
from the stationary state. Very importantly, the pseu- 
doinverse 1Z of the Liouvillean, defined as Q£ _1 Q = 1Z 
is tractable by similar numerical methods as used in 
the evaluation of the current (we use the general- 
ized minimum residual method (GMRes) ^I]). Before 
showing results for the current and noise, we discuss 
an important visualization tool. 

4. Wigner functions 

We have found that Wigner functions are an excellent 
interpretative tool for the numerical results obtained 
for the stationary density matrix. The intuitive pic- 
ture comes from the well-known results in the classical 
limit: the Wigner representation (or, equivalently, the 
phase-space representation) of a regularly moving har- 
monic oscillator is a circle. On the other hand, irregu- 
lar motion under the influence of external noise gives 
rise to a Gaussian probability distribution centered at 



the origin. Since the QD can be either empty or oc- 
cupied, it is advantageous to introduce charge resolved 
Wigner functions (n = corresponds to an empty dot, 
while n = 1 represents the occupied dot), defined as 

w n jx, n-fclk <*-§ wr !*+!> 

(7) 

An example of the empty dot Wigner function is given 
in Fig. The shape is consistent with physical intu- 
ition. Consider, for example, positive x and positive p 
(the shuttle is approaching the drain contact). Then 
the probability of an empty dot is large, because the 
extra charge is very likely to leave the shuttle as the 
drain contact is approached. Analogously, at negative 
x and positive p the probability of an empty dot is 
very small, because the dot has been recharged in the 
vicinity of the left (source) contact, and the charge 
cannot have left the QD because of the exponentially 
small coupling to the right (drain) contact. Very inter- 
estingly, in Fig. 2 we also see an enhanced probability 
located at the origin of the phase space. The inter- 
pretation is that at these values of parameters there 
are two co-existing transport regimes: (i) the charge 
shuttling regime (represented by the ring), and (ii) an 
incoherent tunnelling regime in which charges tunnel 
into and out from the QD uncorrelated to its position. 




Figure 2: The Wigner representation of the empty-dot den- 
sity matrix in the co-existence regime. 



5. Numerical results 

Figure 3 shows the computed stationary current as a 
function of the damping rate. Noteworthy features are: 
(i) One observes a shuttling transition (or, perhaps 
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more appropriately, a cross-over) from a low value of 
current (tunnelling regime) to a high value of current 
(shuttling regime) even in the quantum regime (small 
A); (ii) in the low damping limit the current saturates 
to a universal value of 1 /2tt corresponding to precisely 
one electron transmitted per cycle. 

Current 
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Figure 3: Stationary current as a function of damping rate. 

More dramatic effects are observed in the Fano factor, 
given in Fig. 4. Again, we just list the essential fea- 
tures: (i) For high values of damping, the noise has its 
Poissonian value (~ 1); (ii) At low values of damping, 
the noise is very low (reflecting the orderly nature of 
the shuttling regime); (iii) At the shuttling cross-over 
there occurs a large enhancement of the noise. 




y/co 

Figure 4: The Fano factor as a function of the damping. 



The extraordinarily large values of the Fano factor of 
the order of 600 can be explained as being a conse- 
quence of a slow switching process between two com- 
peting current channels (shuttling and tunnelling), 
and in Ref. 6 we give a detailed analysis of this 
phenomenon, also supported by semianalytic consid- 
erations. 

In conclusion, we have presented a numerical tech- 
nique for solving the generalized Master equation gov- 
erning a generic quantum nanoelectromechanical de- 
vice. The obtained numerical results are interpreted 
with the help of phase space representations. We be- 
lieve that the methods discussed here are also applica- 
ble to many other quantum transport situations, where 
the matrix representations of the relevant operators 
are very large, but where only certain extremal eigen- 
values are important. 

References 

Elements of this work have been summarized in an extended 
abstract published by IEEE. 

1. H. Park, J. Park, A. K. L. Lim, E. H. Anderson, A. P. 
Alivisatos, and P. L. McEuen, Nature 407, 57 (2000). 

2. L. Y. Gorelik, A. Isacsson, M. V. Voinova, B. Kasemo, R. I. 
Shekhter, and M. Jonson, Phys. Rev. Lett. 80, 4526 (1998). 

3. T. Novotny, A. Donarini, and A. P. Jauho, Phys. Rev. Lett. 
90, 256801 (2003). 

4. T. Novotny, A. Donarini, C. Flindt, and A. P. Jauho, Phys. 
Rev. Lett. 92, 248302 (2004). 

5. C. Flindt, T. Novotny, and A. P. Jauho, to appear in Phys. 
Rev. B (2004); cond-mat /0405512| 

6. C. Flindt, T. Novotny, and A. P. Jauho, to appear in Euro- 
phys. Lett. (2004); cond-mat/0410322| 

7. A. D. Armour and A. MacKinnon, Phys. Rev. B 66, 035333 
(2002). 

8. We show in |3 El El that the off-diagonal elements Poi/io 
decouple from Poo/n? an d are not needed in the evaluation 
of the stationary limit. 

9. A good introduction to the Arnoldi scheme can be found 
in G. H. Golub and C.F. Loan, Matrix Computations, The 
Johns Hopkins University Press, 3rd Edition (1996). 

10. Prof. T. Eirola, Private communication. 

11. C. Flindt, Ma ster's Thesis, MIC, Technical Un iversity of 
Denmark, URL http://www.mic.dtu.dk/research /Theoret- 
icalNano /publications/ Theses.htm 



